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ABSTRACT 


This report investigates the utility of the Hilbert-Huang transfonn for the analysis of 
aeroelastic flight data. It is well known that the classical Hilbert transform can be used for 
time-frequency analysis of functions or signals. Unfortunately, the Hilbert transform can only 
be effectively applied to an extremely small class of signals, namely those that are characterized 
by a single frequency component at any instant in time. The recently-developed Hilbert-Huang 
algorithm addresses the limitations of the classical Hilbert transform through a process kn own as 
empirical mode decomposition. Using this approach, the data is filtered into a series of intrinsic 
mode functions, each of which admits a well-behaved Hilbert transform. In this manner, the 
Hilbert-Huang algorithm affords time-frequency analysis of a large class of signals. This powerful 
tool has been applied in the analysis of scientific data, structural system identification, mechanical 
system fault detection, and even image processing. The purpose of this report is to demonstrate the 
potential applications of the Hilbert-Huang algorithm for the analysis of aeroelastic systems, with 
improvements such as localized-online processing. Applications for correlations between system 
input and output, and amongst output sensors, are discussed to characterize the time-varying 
amplitude and frequency correlations present in the various components of multiple data channels. 
Online stability analyses and modal identification are also presented. Examples are given using 
aeroelastic test data from the F- 1 8 Active Aeroelastic Wing airplane, an Aerostructures Test Wing, and 
pitch-plunge simulation. 


NOMENCLATURE 


a 

constant 

a j 

Fourier coefficient 

A x (t),Ay(t) 

time- varying amplitude of Z x (t),Zy(t) 

ATW 

Aerostructures Test Wing 

b 

chord length 

Ci(t) 

intrinsic mode function 

a 

plunge and pitch damping 

Cl a ,C,n a 

lift and moment coefficients due to pitch 

C ' l P ’ Cm (i 

lift and moment coefficients due to control surface rotation, /i 

ELCC 

empirical local correlation coefficient 

EMD 

empirical mode decomposition 

FM 

frequency modulation 

h 

plunge degree-of-freedom 

h(t) 

impulse response function 

h(t) 

Hilbert transfonn of the impulse response function 

h(w) 

marginal spectrum 


HHT 


Hilbert-Huang transform 

HLCC Hilbert local correlation coefficient 

H(co,t) Hilbert spectrum 

2 

H ( (0, t) energy spectrum 

H x (ri x ,t),H y (Vy>t) Hilbert empirigrams from Z x (t),Z y (t) 

H xy (rj,t) Hilbert cross-empirigram of H x (Tj,t),H y (Tj,t) 

IE instantaneous energy 

IM instantaneous magnitude 

IMF intrinsic mode function 

IP instantaneous phase 

ITF instantaneous transfer function 

I a pitch moment-of-inertia 

k h , k a plunge and pitch stiffness 

m mass 

M complex matrix operator of input analytic IMF’s to output analytic IMF’s 

01 orthogonality index 

PV Cauchy principal value 

r(t) residual of empirical mode decomposition 

STD common standard deviation 

SVD singular value decomposition 

U flow velocity 

x(t) real-time domain signal 

x H (t),y H (t) Hilbert transform of x(t),y(t) 

x elastic axis location 

Z x (t),Z y (t) complex analytic signal from Hilbert transform of x(t),y(t) 

Z xy (t) cross-analytic signal from Z x (t),Z y (t) 

a pitch degree-of-freedom 

1 3 control surface rotation 

C damping ratio 

77 intrinsic mode function number 

f] x ,r] x intrinsic mode function numbers from signals x,y 

jjA(M) structured singular value with uncertainty block structure A 

p air density 
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p[ ] 

spectral radius (largest magnitude eigenvalue) 

°xy 

77 -vector of singular values from the empirigrams H x (ij,t),Hy( i]j ) 

a 

maximum singular value 

0 x (t),0y(t) 

time-varying phase of Z x (t),Z y (t) 

Qox (Oi $0y (t) 

initial phase of Z x (t),Z y (t) 

0xy (0 

cross-analytic phase of Z^t) 

CO(t) 

instantaneous frequency 


damped natural frequency 

(O' 

Fourier discrete frequency 


natural frequency 


INTRODUCTION 

The Hilbert transform is a classical tool that has been used in the structural dynamics community 
as an indicator of nonlinearity. It has also been used to estimate nonlinear damping and stiffness 
functions for single degree-of-freedom systems. The Hilbert transfonn expresses a signal as a 
harmonic with time-dependent frequency and amplitude. In this respect, it is an ideal tool for the 
analysis of nonstationary data. Unfortunately, the Hilbert transform has several shortcomings that 
limit its usefulness in practice. Most notably, the Hilbert transfonn computes a single instantaneous 
frequency for a signal at each instant in time. Therefore, when applied to a multi-component 
signal (that is, a signal from a system with multiple modes), the Hilbert transfonn computes an 
instantaneous frequency that conesponds to a weighted average of the component frequencies. 
Such an instantaneous frequency does not provide any information as to the values of the individual 
component frequencies (ref. 1). A further limitation is that the Hilbert transform yields grossly 
distorted estimates of the frequency when applied to signals with nonzero mean and signals which 
have more extrema than zero crossings. 

In order to address these shortcomings, an empirical mode decomposition (EMD) was developed 
by Huang et al. (refs. 2 and 3) as a means of decomposing a signal into a series of components known 
as intrinsic mode functions (IMFs). These IMFs are computed based on local characteristics of the 
signal and can be viewed as an adaptive, data-dependent basis. The IMFs form a complete, nearly 
orthogonal set of basis functions. Most importantly, each IMF contains only a single frequency 
component at any instant in time and therefore admits a well-behaved Hilbert transfonn. Taken 
collectively, the Hilbert spectra of the IMFs yield complete time-frequency information about the 
original signal. This approach, which has been termed the Hilbert-Huang Transform (HHT), makes 
it possible to apply the Hilbert transform to an extremely general class of functions and signals. 

Although a relatively new tool, the HHT algorithm has received considerable attention in a 
number of engineering disciplines. This HHT algorithm has been applied in the analysis of scientific 
data (refs. 2-4), structural system identification (refs. 5-7), and mechanical system fault detection 
(refs. 8-11). A recent adoption into the image processing field (refs. 12-14), the two-dimensional 
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EMD is an adaptive image decomposition without the limitations from filter kernels or cost 
functions. The IMFs are interpreted as spatial frequency subbands with varying center frequency 
and bandwidth along the image. The EMD is a truly empirical method, not based on the Fourier 
frequency approach but related to the locations of extrema points and zero crossings. Based on this, 
the concept of “empiquency,” used for time or space and short for “empirical mode frequency,” 
was adopted to describe signal (image) oscillations based on the reciprocal distance between 
two consecutive extrema points. Fligh concentrations of extrema points have high empiquency 
with sparse areas having low empiquency. Hence, applications for time-frequency-space signal 
processing are feasible. 

A problem with the very versatile and most commonly used Morlet wavelet in dynamics data 
analysis is its leakage generated by the limited length of the basic wavelet function, which makes 
the quantitative definition of the energy-frequency-time distribution difficult. The interpretation of 
the wavelet can also be counterintuitive. For example, definition of a local event in any frequency 
range requires analysis in the high-frequency range, for the higher the frequency the more localized 
the basic wavelet. A local event in the low-frequency range requires an extended period of time 
to discern it. Such interpretation can be difficult if possible at all. Another problem with wavelet 
analysis is its nonadaptive nature. Wavelet basis functions are predefined, whether one type, or a 
multi-wavelet basis, or a dictionary of wavelets is selected. Data analysis is then constrained to 
these bases. Since the Morlet wavelet is Fourier based, it also suffers the many shortcomings of 
Fourier spectral analysis. A truly adaptive basis is a necessary requirement for nonstationary and 
nonlinear time series analysis and should be based on and derived from the data. 

In this report, results are obtained from the new approach. With the HHT, the IMFs yield 
instantaneous frequencies as functions of time that give sharp identification of imbedded structures. 
The main conceptual innovation in this approach is the introduction of the instantaneous frequencies 
for complicated data sets, which eliminate the need of spurious harmonics to represent nonlinear and 
nonstationary signals. This report looks at the effect of enhancements like local-on-line versions 
of the algorithm (ref. 15). To date, HHT analysis has only been performed on individual signals 
without regard to correlation with other data channels or system inputs-to-outputs. Application 
for correlations between system signals are introduced to characterize the time-varying amplitude 
and frequency modulations present in the various components of multiple data channels including 
input and distributed sensors. In these respects, this report attempts to elucidate the way EMD 
behaves in the analysis of the F-18 Active Aeroelastic Wing (AAW) airpane (ref. 16), aeroelastic 
and aeroservoelastic flight test data, as well as Aerostructures Test Wing (ATW) (ref. 17) and 
pitch-plunge simulation data. 

EMPIRICAL MODE DECOMPOSITION 

The classical Hilbert transform is, in principle, an effective tool for time-frequency analysis. 
Unfortunately, in practice, it can only be applied to an extremely restricted class of signals. In 
order to extend the utility of the Hilbert transform to more general signals, Huang et al. (refs. 2 
and 3) developed the EMD as a means of preprocessing data before applying the Hilbert transform. 
The EMD procedure decomposes the original signal into a set of IMFs, each of which admits a 
well-behaved Hilbert transfonn. The Hilbert transform is then applied to each individual IMF, 
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yielding an instantaneous frequency and amplitude for each IMF. Therefore, this procedure, 
termed the Hilbert-Huang transform (HHT), enables the Hilbert transform to be applied to 
multi-component signals. 

There are two criteria that each IMF must satisfy in order to be amenable to the Hilbert transform. 
Namely, each IMF must have zero mean and the number of local extrema and zero crossings in 
each IMF can differ by no more than one. The first step in the EMD procedure is to connect all 
the local maxima of the original signal using a cubic spline. Similarly, the local minima are also 
connected with a cubic spline. The two splines define the envelope, and the mean of this envelope 
is then calculated and subtracted from the original signal. The resulting signal is then tested to see 
if it satisfies the criteria for an IMF. If it does not, the sifting process is repeated until a suitable 
first IMF, ci(t), is obtained. The sifting process is then applied to the residual signal x(t)~ q (/) 
to obtain the next IMF. This process is repeated until all that is left is a final residual, r(t ) , which 
represents the trend in the data and is not an IMF. Therefore, as shown in equation (1), the EMD 
procedure yields a decomposition of the original signal in terms of n IMFs and the residual. 


x(t) ='£cj(t) + r(t) (1) 

7=1 


The number of IMFs obtained is dependent on the original signal. 

The EMD procedure serves to generate IMFs that are amenable to the Hilbert transform. In 
particular, the manner in which the EMD is perfonned guarantees that each IMF will only possess 
a single harmonic at any instant in time. Currently, the EMD procedure is ad hoc in the sense that 
there is no rigorous mathematical theory behind it. Recent attempts have been made to formalize 
EMD, most notably the work of those who explored the properties of B-splines and their use in the 
EMD process. Despite the lack of a firm theoretical foundation, several mathematical properties of 
the IMFs are well understood. In particular, it is clear from equation (1) that the IMFs, along with 
the residual, fonn a complete basis for the original signal x(t ) . In addition, Huang et al. (ref. 2) 
demonstrated that the IMFs are nearly orthogonal with an orthogonality index [OI m „ defined in 
equation (2)] for the IMFs. This definition seems to be global but actually only applies locally. 
Adjacent IMFs could have data with the same frequency but at different times. 
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mn 


m n c (t)c k (t) 

Z-i Z-i Z-i . . 2 

t=0 j=\k(*j )= 1 


( 2 ) 


The EMD does not yield a unique basis for the original signal since there are countless sets 
of suitable IMFs that can be generated from a given signal. However, the various IMF sets from 
the different sifting criteria are all equally valid representations of the data provided their OIs are 
sufficiently small (refs. 2 and 4). The IMFs depend on the stoppage criterion, maximum number 
of siftings, intermittance criteria, the end point boundary conditions, and use of curvature- or 
extrema-based sifting (refs. 2-4). The uniqueness problem can only be meaningful if all these 
parameters are fixed a priori. The problem is how to optimize the sifting procedure to produce 
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the “best” IMF set. These questions are difficult to answer theoretically. In Huang et al. (ref. 4), a 
confidence limit is defined for the first time without invoking the ergodic assumption. This provides 
a stable range of stopping criteria for the EMD-sifting operation, thereby making the HHT method 
more definitive. Discussion later points out that the EMD sifting process also acts in such a manner 
as to obtain IMFs that correspond to approximate bandpass filtering of the original signal. 

Once the EMD procedure has been used to generate a set of IMFs, the Hilbert transform can 
be applied to each individual IMF. In this manner, an instantaneous frequency and amplitude is 
computed for each function. A common method for displaying the Hilbert spectrum is to generate 
a two-dimensional plot with time and frequency axes. The amplitude is then plotted as a color 
spectrogram in the time-frequency plane. By plotting the Hilbert spectra of all the IMFs together 
one obtains complete time-frequency information about the original signal. 

To demonstrate the importance of the EMD in obtaining meaningful Hilbert spectra for general 
signals, consider the following signals: 

X] (t) = sin( 2 nt) + sin( 8 nt) + sin( 32m) 
x 2 (t) = sin(4m) + t 

The signal Xj is a combination of harmonics with frequencies of {1, 4, 16} Hz while the signal 
x, is composed of a ramp signal and a 2-Hz sine wave. First, consider the Hilbert spectra obtained 
by directly applying the Hilbert transform to each signal. Figure 1 displays the resulting Hilbert 
spectra for both signals. The Hilbert spectrum of x L yields a fluctuating frequency with most of 
the energy concentrated at 7 Hz. This corresponds to the average of the three frequencies present 
in Xj (each of which has the same amplitude) and gives no useful information about the original 
signal. This averaging is due to the fact that the Hilbert transform computes a single instantaneous 
frequency for the signal. The Hilbert spectrum of x, illustrates the effect of a nonzero mean. In this 
case, the estimated frequency is distorted and the 2-Hz harmonic is not identified at all. 

Figure 2 displays the signal x { and its decomposition into four IMFs (imf 1 = c { , imf 2 = c 2 , etc.) 
and a residual. Clearly, each IMF has captured a different frequency component of the original 
signal. The fourth IMF is extremely small and results from boundary effects in the EMD process. 
The Hilbert spectra of the individual IMFs are plotted together in figure 2. Now, an effective 
Hilbert spectrum, or Hilbert-Huang spectrum, of x l has been obtained, with all three frequency 
components clearly identified. Similarly, figure 3 depicts two IMFs and a residual generated 
from the signal x 2 . The 2-Hz frequency component has been sifted into the first IMF and the 
ramp component has been identified as the second IMF. The mean of the signal is 2, which has 
been separated out as the residual. The Hilbert spectrum clearly shows the 2-Hz component and 
estimates a low frequency for the ramp component. This occurs because the ramp is treated as part 
of a low frequency wave. The IMFs and Hilbert spectra in figures 2 and 3 illustrate that there are 
some minor boundary effects associated with the EMD process. Most importantly, these examples 
demonstrate that EMD makes it possible to apply the Hilbert transform to signals that otherwise do 
not admit well-behaved Hilbert transforms. 
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Finally, two more examples (ref. 15) illustrate automatic and adaptive (signal-dependent) 
time-variant filtering of general mixtures of signals. A signal composed of three components which 
significantly overlap in time and frequency is successfully decomposed in figure 4 for the sum of 
two sinusoidal frequency modulations and one Gaussian wave packet. Another example, accenting 
the nonharmonic nature of EMD, is given in figure 5. The analyzed signal (top) is the sum of 
three components, a sinusoid superimposed on two triangular waveforms with periods smaller 
and larger than the sinusoid. The decomposition performed by the EMD is given in the three 
IMFs and the residue. In this case, both linear (sinusoid) and nonlinear (triangular) oscillations are 
effectively identified and separated. Flarmonic analysis (Fourier, Morlet wavelets) would produce 
a less succinct and more nondescriptive decomposition. 

Filtering Properties 

The filtering properties of EMD have been studied in some detail (refs. 18-20). The EMD 
process yields a data-dependent decomposition that focuses on local characteristics of the signal. 
In particular, EMD sifts out the highest-frequency component in the signal at any given time. 
Indeed, EMD has surprisingly been shown to behave as a dyadic bandpass filter when decomposing 
Gaussian white noise, much like a multiresolution wavelet decomposition. However, the EMD 
method as an equivalent dyadic filter bank is only in the sense of its global behavior over the entire 
time extent. In representing the time-frequency distribution, the Hilbert spectrum of each IMF is 
actually localized at any time. This is different from predetermined filtering such as with Morlet 
wavelet filtering. 


Analytical Interpretation 

The EMD is faced with the fundamental difficulty of not admitting an analytical definition, but 
of rather being defined by the algorithm itself, thus making the analysis of its performance and 
limitations difficult. The need for rigorous mathematical foundation is imperative. This fundamental 
problem of the EMD has to be resolved since only with the IMF can nonlinear distorted waveforms 
be resolved from nonlinear processes. There have been attempts to circumvent the mathematical 
difficulties in the EMD with some success by casting the IMFs in terms of B-splines (ref. 21), 
and applying toward mechanical system fault-detection (ref. 9). System identification of the IMFs 
as a multi-component system is suggestive in the light of multiresolution system identification 
procedures such as multiresolution singular value decompositions (SVDs), Kalman filters, and 
subspace algorithms. 

More significantly, characterization of IMFs as solutions to certain self-adjoint ordinary 
differential equations is demonstrated (refs. 22 and 23). Construction of envelopes which do not rely 
on the Hilbert transform is used directly to compute the coefficients of the differential equations. 
These equations are natural models for linear vibrational problems and provide further insight 
into both the EMD procedure and utilizing its IMF components to identify systems of differential 
equations naturally associated with the components. One of the uses of the EMD procedure is 
to study solutions to differential equations, and vibration analysis was a major motivation in the 
development of the Sturm-Liouville theory. 
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HILBERT TRANSFORM AND INSTANTANEOUS FREQUENCY 


The Hilbert transform of a time-domain function or signal x is defined in equation (3), 

] r°° x(T^) 

y(t) = H{x(t)} = -PV\ — —dz (3) 

where PV denotes the Cauchy principal value, needed because the integrand is singular at z = 1 . 

The Hilbert transform can be viewed as the convolution of the original signal with 1 It, 
emphasizing temporal locality of x(t) . Note that, unlike Fourier analysis, the Hilbert transform of 
a time-domain signal is another time-domain signal. The Hilbert transfonn is sometimes applied 
to frequency-domain signals using a similar expression as equation (3), but this report will focus 
on the time-domain case. In practice, the Hilbert transform is usually calculated using the Fourier 
transfonn. Therefore, the fast Fourier transfonn algorithm can be employed for the efficient 
calculation of the Hilbert transform. 

A signal, x, and its Hilbert transform, y, can be used to define a complex analytical signal as in 
equation (4). 


z(t ) = a(t )e l ^ ( ' t ^ = x{t) + iy(t ) 


( 4 ) 


Therefore, the Hilbert transform pair \x(t),y(t ) } can be expressed as a harmonic function with 
time -varying amplitude a(t) and time-varying phase angle 0(t). 


ait ) = y]x(t) 2 +y(tY 


Q(t) = tan 


-1 


y( 0 
40. 


Given the time-dependent phase angle, the instantaneous frequency of the signal can be defined 
as (ref. 24) 


co(t) = 


ddjt) 

dt 


In the context of equation (3) with instantaneous frequency, the Hilbert transfonn of an 
IMF can be interpreted as giving the best fit with a sinusoidal function to the data weighted by 
1 It. Instantaneous frequency can be computed using the derivative definition or centralized 
finite difference (ref. 24). In general, there are an infinite number of ways to express a signal 
as in equation (4), so there can also be an infinite number of instantaneous frequencies. The 
Hilbert-transfonn pair was proposed to uniquely define the amplitude and phase by building the 
complex analytic signal from the given signal with the original signal, x(t), as the real part and the 
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orthogonal transformed signal, y(t), the imaginary part, out of phase with x(t) by 
Fourier analysis yields the decomposition 


— . Note that 
2 


oo 



which is similar to the form of the Hilbert transform in equation (4). A key difference, however, is 
that the Fourier decomposition is in terms of harmonics with constant amplitudes and frequencies. 
In contrast, the Hilbert transfonn yields instantaneous amplitudes and frequencies. Therefore, in 
principle, the Hilbert transform is an ideal tool for the time-frequency analysis of a general class of 
signals, including nonstationary signals. 

Another important distinction between Fourier analysis and the Hilbert transform is that the 
Fourier decomposition is in terms of multiple harmonics of constant amplitude and frequency, 
thereby producing artificial harmonics to maintain energy conservation for nonstationary and 
nonlinear data. In contrast, the Hilbert transfonn of a signal yields an expression in terms of a single 
harmonic with a time-varying frequency and amplitude. For this reason, the Hilbert transform 
is only suitable for the analysis of mono-component signals, or signals that are composed of a 
single frequency component at any instant in time. This is a considerable limitation as it implies 
that the Hilbert transform cannot be directly applied to signals that are composed of multiple 
harmonics. As was shown, the Hilbert transform fails to identify the individual frequencies and 
instead computes a single instantaneous frequency that corresponds to a weighted average of 
the component frequencies. The resulting instantaneous frequency is both physically invalid and 
erroneous (ref. 1) since a multi-component signal has more than one instantaneous frequency. An 
additional limitation of the Hilbert transform is that it yields extremely distorted estimates of the 
instantaneous frequencies of signals with nonzero mean and signals that have more local extrema 
than zero crossings. 

The EMD responds to the dilemma surrounding the applicability of instantaneous frequency. It 
decomposes a multi-component signal into its associated mono-components while not obscuring or 
obliterating the physical essentials of the signal and allows the traditional definition of instantaneous 
frequency to be complete by being applicable to signals of both mono- and multi-component. 
To follow the true frequency evolution within a multi-component signal, it is necessary to break 
down the components into individual and physically meaningful intrinsic parts. The adaptive 
and nonarbitrary decomposition using EMD produces an orthogonal set of intrinsic components 
each retaining the true physical characteristics of the original signal. The mono-components or 
intrinsic modes satisfy the conditions for a well-defined notion of instantaneous frequency. These 
conditions include symmetry, no dependence on predefined time scales, revelation of the nature of 
simultaneous amplitude and phase variation, and near-orthogonality. 

Finally, since instantaneous frequency displays frequency variation with time, changes of 
dynamic states indicative of nonlinearity can be identified. For example, if the instantaneous 
frequency of a new mode is about half of the frequency of the old mode in a bifurcation, period 
doubling occurs. If the instantaneous frequency of the new mode is disproportionate with the 
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old mode, quasi-periodic bifurcation occurs. Similarly, intermittence and chaotic motion can be 
determined. A dynamic state can be diagnosed simultaneously by observing the changes in time 
of the instantaneous frequency components and their corresponding energy. In summary, the 
concept of mode defined in relation to instantaneous frequency as a periodic-modal structure in 
the instantaneous time-frequency plane is found to be more appropriate than artificial sinusoidal 
harmonics in characterizing nonlinear responses (ref. 10). Instantaneous frequency is a quantity 
critical for understanding nonstationary and nonlinear processes. 

LOCAL-ON-LINE DECOMPOSITIONS 

In the original EMD formulation, sifting iterations are applied to the predefined full length 
signal as long as there exists a locality at which the mean of the upper and lower envelopes is not 
considered sufficiently small enough. Excessive iterations on the entire signal to achieve a better 
local approximation contaminates other parts of the signal by overcompressing the amplitude 
and overdecomposing it by spreading out its components over adjacent intrinsic modes, that is, 
overiteration leads to overdecomposition. The various stoppage criteria (to fulfill that the number 
of extrema and the number of zero crossings must differ at most by one, and that the mean 
between the upper and lower envelopes must be close to zero) are attempts to avoid the rigor of 
the symmetry of the envelopes without a mathematically rigorous definition for an adaptive basis. 
The hierarchical and nonlinear nature of the EMD algorithm will not provide that the EMD of 
segmented signals will be the segmentation of individual EMDs. Therefore, a variation referred to 
as local EMD (ref. 15) introduces an intermediate step in the sifting process. Localities at which 
the error remains large are identified and isolated, and extra iterations are applied only to them. 
This is achieved by introducing a weighting function such that maximum weighting is on those 
connected segments above a threshold amplitude, with a soft decay to zero outside those supports, 
much like soft-thresholding is done in wavelet denoising. 

Another option is based on the idea that sifting relies on interpolations between extrema, and thus 
only requires a finite number of them (five minima and five maxima in the case of the recommended 
cubic splines) (refs. 2 and 21) for local interpolation. Extraction of a mode could therefore be 
moving blockwise instead of globally over the entire time span. This led to the development of 
the on-line EMD (ref. 15). A prerequisite for the sliding window extraction of a mode is to apply 
the same number of sifting steps to all blocks in order to prevent possible discontinuities. Since 
this would require the knowledge of the entire signal, the number of sifting operations is proposed 
to be fixed a priori to a number less than 10 for effective application of the on-line version of the 
EMD algorithm operating in coordination with the local EMD described above. The leading edge 
of the window progresses when new data become available, whereas the trailing edge progresses 
by blocks when the stopping criterion is met on a block. Therefore, an IMF and corresponding 
residual are computed sequentially, then again applied to this residual, thus extracting the next 
mode with some delay. 

An aileron command multisine input used on the F-18 AAW for aeroservoelastic response and 
flutter clearance is shown in figure 6 (top plots) using the standard EMD (left) and local-on-line 
version (right). The bandpass nature of IMFs (refs. 18-20) is reflected in the three standard IMF 
mean frequencies, {23.9, 16.3, 7.7}Hz for each of IMFs {1#, #2, #3}, respectively, and on-line 
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corresponding IMF mean frequencies {23.6, 13.0, 7.2} Hz. Immediately noticeable is the more 
efficient extraction of the signal components by the local-on-line algorithm, most evident by the 
sparse second and third IMFs (imfZ and imO) being more sparse than the corresponding standard 
IMFs. Besides the obvious advantage of an on-line algorithm for decomposing data, it has been 
found to clearly surpass the standard global algorithm in terms of computational burden, especially 
with long original data records. An added bonus is that it generally has better orthogonality 
properties among the IMFs, witnessed by an order of magnitude improvement in the 01 defined 
in equation (2). 


LOCAL ANALYTIC SIGNAL CORRELATION 

Since the IMFs allow permissible, meaningful, physically sensible, and unique interpretations 
of instantaneous frequency of general signals of interest, they fit into the class of asymptotic signals 
such that the time variation of the IMF amplitude and frequency may directly be recovered from 
the time variation of the amplitude and of the phase derivative of the associated analytic signal. 
An IMF after performing the Hilbert transform can be written as in equation (4). These complex 
components are now used for analysis of analytic data correlations (refs. 25 and 26) between 
input-output and amongst spatially distributed sensor outputs. Note that these analyses are all 
local in nature since there are no assumptions of stationarity, and differ from classical double-time 
expressions (ref. 25) by being instantaneous. 

Local Correlation Coefficient 

Correlations are made between transformed IMFs of various signals given the associated 
complex analytic signals 


Z x (t) = A x (t)e i ^ (t) =x(t) + ix H (t) 
Z y ( t) = Ay ( t)e l<t>y (t) = y(t) + iy H ( t) 


by considering the cross-analytic signal defined in equation (5) by 


* id) (t) i (^v 00 ) 

Z^(t) = ZJt)Z y (t) = AJt)A y (t)e l ^ (t) =AJt)Ay(t)e L y ° y °* J 


( 5 ) 


which is not necessarily analytic but can be used to characterize the time variation of the phase 
difference between components at similar frequencies extracted from two simultaneously acquired 
signals. This Hilbert demodulation technique is only applicable for analytic components, in this 
case the transformed IMFs of the original signals. If the two components are mainly in-phase and 
a change occurs in the cross-analytic signal, the instantaneous frequencies of the two components 
may be drifting apart and this produces an instantaneous contribution to the phase difference from 
equation (5). A measure of the local correlation between components, in terms of simultaneous 
changes in instantaneous amplitude or frequency (phase) between analytic signals, is the Hilbert 
local correlation coefficient, HLCC. 
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HLCC(t) = 



As a simple illustration of the utility of the HLCC, figure 7 shows the EMDs of wingtip 
accelerometer responses (accels) due to a symmetric aileron input, with the top plot in each set 
being the original signal, followed by EMDs 1 to 3 and the residual. Figures 7(a) and 7(b) represent 
the forward (fwd) wingtip accels, and corresponding IMFs seem not to compare quite as well as 
the IMF sets of aft wingtip accels, shown in figures 7(c) and 7(d), for the same input. The IMF sets 
shown in figures 7(c) and 7(d) seem much more correlated with each other than those shown in 
figures 7(a) and 7(b). 

Table 1 lists the means, medians, and standard deviations of the HLCC data over the time span of 
the maneuver. Note that HLCC indicates correlation coefficients between Accel#m and Accel#n, 
using cross-analytic IMF function Z correlating Accel#m IMFx to Accel#n IMFv. Results from 
the table indicate strongest HLCC s for Z when x = y, as expected since this is a correlation 
between similar bandpass characteristics of IMFs from different responses due to the same input. 
For comparison, the standard statistical signal correlation coefficient between Accel# 1 and Accel#3 
is C 13 = 0.8633, and between Accel#2 and Accel#4 is C 24 = 0.8332 (correlating fwd-to-fwd accels 
and aft-to-aft accels, respectively, wingtip-to-wingtip). This is consistent with the original signal 
plots in figures 7(a) and 7(b) in that the forward wingtip accels seem to correlate reasonably well 
with each other, as also between the aft wingtip accels, but not forward compared with aft. Evident 
in table 1 is the common standard deviations (STD in all the correlations, near STD = 0.7). In the 
HLCC values good commonality is found between the mean and medians comparing HLCC 13 
and HLCC 24 in that the trend is similar amongst the cross-analytic IMF functions Z xy where x 
and y correspond to respective IMF numbers in different accel EMDs. This is surprising, since in 
figure 7 the correlation in corresponding IMFs seems much worse between the two fwd accels than 
between the two aft accels. 


Table 1. Hilbert local correlation coefficient results from F/A-18 active aeroelastic wing wingtip 
acceleration data. 


Z *y 

Mean HLCC 13 

Median HLCC 13 

STD 13 

Mean HLCC 24 

Median HLCC 24 

STD 24 

Z n 

0.3283 

0.5968 

0.6801 

0.2733 

0.5309 

0.7060 

Z 12 

0.0917 

0.1782 

0.7171 

0.0502 

0.1080 

0.7142 

Z 13 

0.0098 

0.0252 

0.7113 

0.0005 

-0.0106 

0.7057 

Z 2 1 

0.0953 

0.1870 

0.7148 

0.0492 

0.0896 

0.7095 

Z 22 

0.1844 

0.3500 

0.7048 

0.2197 

0.4012 

0.7023 

Z 2 3 

0.1187 

0.2364 

0.7136 

0.0692 

0.1319 

0.7052 

Z31 

-0.0071 

-0.0175 

0.7075 

0.0032 

0.0062 

0.7056 

Z 3 2 

0.0845 

0.1753 

0.7097 

0.0996 

0.2025 

0.7183 

Z33 

0.2403 

0.4425 

0.6866 

0.1553 

0.2877 

0.7027 
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Figure 8 shows the HLCC functions between Accel#l and Accel#3 ( HLCC' xy ) at 5 to 20 s. 
Decomposing and representing correlations in this fashion allows a true time-localized instantaneous 
measure revealing subtle properties in the data and inconspicuous relations to other data sets using 
analytic components. 


Instantaneous Transfer Function 

In addition to investigating correlations between sensors, an instantaneous system transfer 
function (ITF) is introduced with instantaneous magnitude ( IM) and the HLCC as a phase 
parameter between input-output analytic signals. This is viable in terms of the HHT inducing 
analytic properties to the IMFs to yield localized system input output properties. Instantaneous 
transfer function, its IM, and its instantaneous phase (IP) are defined. 

Z (t) 

ITF(t) = — — \lM(l) = \lTF(t)\;IP(t) = cos -1 \ HLCC (I)} = (p xy (l) 

Zxx(t) 

Figure 9(a) represents IM (IM c ( y ) and figure 9(b) represents IP (IP xy ) from the ITF ( ITF xy ) of 
the aileron command input (figure 6) to Accel# 1 [figure 7(b)] analytic IMFs (Hilbert-transformed, 
that is, {Z T , Z ( ,}). Note that ITF™ 1 indicates transfer functions between Input#m and Accel#n, 

using cross-analytic IMF function Z vv correlating Input#m’s IMFx to Accel#n’s IMFv. In this 
sense, by interpreting the input-output Hilbert-transformed pairs of IMFs as a multi-component 
system of input-output signals, local stability measures are deemed to be feasible by tracking gain, 
phase, and instantaneous frequencies between each IMF pair. How these local IMF properties 
correspond to global system properties, given the analytic transfonned IMFs, and aeroelastic and 
aeroservoelastic applications in stability and health monitoring, is currently being researched. 

Empirigrams and Empirical Local Correlation Coefficient 

Given two EMDs from different signals, select a common number of IMFs to use for 
correlation (as three were selected in earlier examples for input and outputs). This is generally not 
difficult, especially for aeroelasticity data, since the higher-numbered IMFs approach a residual 
characteristic quickly. From the transformed IMFs,{Z.(t), Z it)}, define the corresponding set of 
two-dimensional Hilbert empirigrams, H x (rj x ,t),H y (ri y ,t), where for each common >1 = = >] } , the 

IMF number from the input EMD corresponds with the output EMD. Define in equation (6) the 
Hilbert cross-empirigram 


H xy (ri,t) = H* x (ri,t)H y (ri,t) 


( 6 ) 


which correlates the respective IMFs from the two Hilbert empirigrams. Empirigrams relate to 
time-scale wavelet scalograms (ref. 24), which relate to time-frequency maps, since scales relate 
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to frequencies in standard wavelet decompositions. Because of the bandpass nature of IMFs 
discussed previously, a similar construction emanates with the HHT. As with wavelet scalograms, 
the real part of the cross-empirigram (co-empirigram) gives the instantaneous contribution of each 
IMF to the correlation between two signals. An empirical local correlation coefficient ( ELCC ), 
ELCC(jj,t ) 5 i s then defined as 


ELCC(r\,t) 


Re 

1 1 

tC! 

i i 

_ 

Re 

1 1 
3 

i i 



H xy (VC) 


H x (TJ,t)\ 

H y (ri,t)\ 


where ELCC{r],t) (between ±1, as the HLCC) gives the instantaneous contribution between 
corresponding IMFs from the two signals to the correlation coefficient. Figure 10(a) shows an 
imaged decomposition plot of the ELCC between the input signal from figure 6 to Accel# 1 in 
figure 7(b), but includes all nine IMFs from the EMDs of the input and output. The first row 
represents the contributions of the first input IMF to the first output IMF, and so on up to the ninth 
IMF. In each IMF row there is much oscillation (higher frequency in the first and lower frequency 
to the ninth) of contributions from corresponding IMFs to the correlation. There are generally 
stronger correlations over longer time spans in the higher-numbered IMFs (lower frequencies), but 
the lower-numbered IMF correlations are less obvious due to the higher frequency content. There 
is a tendency to cycle from high-to-low-to-high [strong(positive)-to-none-to-strong(negative)] 
correlation very rapidly. Figure 10(b) gives a more detailed depiction, between 10 to 11 s, 
demonstrating a rich interplay between correlation of mid-to-lower IMFs (higher frequencies) 
over the shorter time period. Yet another view is presented in figure 1 1 , where the contours are split 
up discretely in three dimensions showing the heavy emphasis in the higher frequencies because 
the contours are more congested in each of the IMF levels. These representations highlight the 
areas of commonality and incongruity between corresponding input-output IMFs. 

From the Hilbert spectrum H(co,t) the energy spectrum H 2 (a),t) gives instantaneous energy 
(IE), 


IE(t) = J co H 2 (co,t)d(0 

which is an indication of the energy fluctuations with time being weighted by the Hilbert spectrum 
localized energy over the entire set of IMFs. Corresponding to the signals in figure 7 are the IE 
profiles of the output accelerometer responses in figure 12(a). The aft wingtip accels (bottom 
plots) are very similar (from the symmetric aileron input), indicating energy at {7-10, 13-20, 30} 
s time locations. In this case, with the particular multisine input from figure 6 programmed over 
the 3 to 35 Hz frequency range, these times correspond closely to the primary F-18 AAW modal 
frequencies in a corresponding frequency range, that is, {6-9, 12-20, 30-35} Hz, as will be shown 
with marginal Hilbert spectra. The forward accelerometer responses in the top plots are also very 
similar indicating modes near {6, 12-17, 25-32} Hz. Instantaneous energy over a sensor suite is 
therefore a nonstationary indicator of time-varying energy distribution amongst the sensor array. 
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MARGINAL SPECTRA 


From the Hilbert spectrum H(a>,t ) also is defined the marginal spectrum. 

h(a>) = ^ t H(a>,t)dt 

Corresponding to the signals in figure 7 are the marginal spectra of the output accelerometer 
responses in figure 12(b). The aft wingtip accels (bottom plots) are again similar showing modal 
response at {5-10, 12-20, 30-35} Hz. The forward accelerometer responses in the top plots are 
also again very similar indicating modes near {6, 12-17, 25-32} Hz. 

The marginal spectrum measures total amplitude (or energy) contributions from each frequency 
value over the entire data record in a probabilistic sense. Frequency in either H(a>,t ) or h(yv) is 
very different from Fourier spectral analysis. In the Fourier analysis the existence of energy at a 
frequency implies a wave component persisting through the data. Energy at the marginal Hilbert 
frequency, however, implies a higher likelihood (expected value over time) for such a wave to 
have appeared locally. As stated earlier, the Hilbert transform of an IMF gives the best fit with a 
sinusoidal function to the data weighted by 1 It, which makes it instantaneous. In fact, the Hilbert 
spectrum is a weighted nonnormalized joint amplitude-frequency-time distribution in which the 
weight assigned to each time-frequency cell is the local amplitude. Consequently, frequency in the 
marginal spectrum indicates only the likelihood that an oscillation with such a frequency exists. 
Exact occurrence of frequency content is given in the full Hilbert spectrum. The Fourier spectrum 
is meaningless for nonstationary data, and there is little similarity between Hilbert and Fourier 
from previous studies (refs. 2 and 3) for nonstationary data. Also, marginal cross-spectra between 
signals does not make sense since the time-dependence (causality) is lost and the frequencies are 
time-ignorant, so an input-output correlation at a certain frequency is meaningless. 

The time-frequency Fourier spectrum, a spectrogram, suffers from the same restrictions over 
time due to windowing distortions. There is a lower bound on the local time-frequency resolution 
uncertainty product of the spectrogram due to the windowing operation. This limitation is an 
inherent property of the spectrogram and is not a property of the signal or a fundamental limit 
(ref. 27). For many other time-frequency distributions, the local uncertainty product is less than 
that of the Fourier spectrogram and can be arbitrarily small. These results are contrary to the 
common notion that the uncertainty principle limits local quantities. Similar considerations apply 
to window-based filter bank (wavelet) methods. This limitation is due to the window and not to any 
inherent property of the signal. 

The fact that windowing is not a factor, and therefore standard time-frequency resolution 
limitations based on the uncertainty principle do not apply, is a key point for application and 
understanding of HHT analysis. As seen in the on-line analysis, an adaptive sliding window method 
has good performance simply requiring adequate data in the window to initiate a sifting process for 
satisfying the two IMF properties for analyticity. 
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WING DYNAMICS ANALYSIS 


Dynamics between different sensors from the same input will now be investigated. This information 
will be used to guide the analysis of a group of wing accelerometer responses from a single input 
in terms of input-output correlation contribution for the F-18 AAW airplane (ref. 16) and ATW 
(ref. 17). 


F-18 Active Aeroelastic Wing Airplane 

Collective F-18 AAW aileron position, used as the multisine input, was obtained as the average of 
four position transducer measurements from the right and left ailerons. Outputs are twelve wing 
structural accelerometers located at the left (six accels) and corresponding right (six accels) outer 
wings, all sampled at 400 sps, as designated in table 2. 

Table 2. F-18 active aeroelastic wing accelerometer nomenclature. 


3. L-wt-fwd 
3. Left-wingtip forward 
7. L-wof-fwd 

7. Left-wing outer-fold forward 
11. L-wif-fwd 

11. Left-wing inner-fold forward 

1. R-wt-fwd 

1 . Right-wingtip forward 
5. R-wof-fwd 

5. Right-wing outer- fold forward 
9. R-wif-fwd 

9. Right-wing inner-fold forward 

4. L-wt-aft 

2. R-wt-aft 

4. Left-wingtip aft 

2. Right-wingtip aft 

8. L-wof-aft 

6. R-wof-aft 

8. Left-wing outer- fold aft 

6. Right-wing outer-fold aft 

12. L-wif-aft 

10. R-wif-aft 

12. Left-wing inner-fold aft 

10. Right-wing inner-fold aft 


Empirical mode decompositions were calculated for the input and all outputs, with mean 
frequencies calculated for all the output IMFs compared with the single input IMFs. Table 3 is a 
summary of the results. Included are the first six IMFs averaged amongst the accels, percentage 
differences of these compared with the input IMFs, and standard deviations of output sensor IMFs. 
The accels correspond very well amongst each other in IMF frequencies (as frequency decreases, 
from left to right). The comparison to the input frequency is excellent even though input and output 
EMDs are performed independently. Orthogonality of the IMFs in each case was also very good. 
This shows that for this type of data, where responses are all from a common input (whether it is 
kn own or not), IMFs can be used for sensor-to-sensor correlation and input-output analysis. 
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Table 3. Intrinsic mode function frequencies from F-18 active aeroelastic wing input collective 
aileron position and outer-wing response data. 


Signal 

IMF1 

IMF2 

IMF3 

IMF4 

IMF5 

IMF6 

0. Input 

116.62 

59.61 

29.34 

16.75 

9.05 

4.90 

1. R-wt-fwd 

99.05 

56.94 

27.90 

12.56 

6.96 

3.62 

2. R-wt-aft 

123.85 

71.36 

35.15 

16.83 

9.35 

5.02 

3. L-wt-fwd 

114.47 

66.05 

33.78 

15.06 

7.95 

5.32 

4. L-wt-aft 

127.56 

78.71 

41.73 

21.06 

11.92 

6.58 

5. R-wof-fwd 

105.83 

58.52 

32.06 

17.17 

9.34 

5.47 

6. R-wof-aft 

96.31 

54.96 

29.24 

14.91 

8.64 

4.63 

7. L-wof-fwd 

115.93 

58.77 

30.88 

16.59 

8.66 

4.32 

8. L-wof-aft 

102.99 

52.76 

27.14 

15.35 

8.65 

4.83 

9. R-wif-fwd 

112.19 

60.61 

31.99 

17.78 

9.90 

5.17 

10. R-wif-aft 

115.39 

64.07 

36.09 

19.05 

10.29 

5.72 

11. L-wif-fwd 

124.91 

66.06 

32.99 

18.32 

10.53 

5.19 

12. L-wif-aft 

129.41 

66.05 

34.02 

17.97 

9.66 

4.92 


Average Output IMF Frequencies 


Absolute 

114.0 

62.9 

32.8 

16.9 

9.3 

5.1 

Wrt Input 

2.3% 

5.2% 

10.4% 

0.8% 

2.9% 

3.2% 


Output IMF Frequency STDs 



Absolute 

11.2 

7.3 

4.0 

2.2 

1.3 

0.7 

Wrt Mean 

9.8% 

11.7% 

12.1% 

13.2% 

13.9% 

14.4% 
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Now the concept of using the essentially-orthogonal analytic IMFs from the inputs and outputs is 
pursued to establish a multi-loop connotation of input IMFs to output IMFs for correlation and even 
stability properties. The input IMFs are interpreted as an orthogonal decomposition of the input(s), 
and the same for output IMFs for output(s). This can be generalized to multi-input-multi-output 
(MIMO) signal analysis where in reality each signal is represented by its EMD. Recall from 
the transformed input-output IMFs, {Z (t), Z(t)}, the corresponding set of empirigrams 
{HJj] x ,t), H y (rj y (t)} were defined in equation (6) for each set o 1' rj = rj x = 77 common IMFs from 
the input and output EMDs. For the analysis described here it is not necessary to use an identical 
number, rj, of IMFs from input and output (but this restriction is maintained here for simplicity). 
Then the Hilbert cross-empirigram H (rj, t ) = H* x (77, t) correlates the respective analytic IMFs 
from the two Hilbert empirigrams. Now define 

0 xy(ri) - cr[//*(77,t)F/J(?7,E)] 

where o ' represents the ^-vector of singular values from the S VD of the product of the empirigrams . 
The singular values represent relative contributions from the principal cross-correlation analytic 
IMFs as a result of correlation of all input analytic IMFs to all output analytic IMFs over the entire 
time span. Therefore, higher a -valued cross-analytic IMFs have more input-output significance 
in terms of operator norm from input to output. The maximum singular value, o = max j; (0 (>])) , 
of this input-output operator corresponds to the structured singular value with a full-complex 
uncertainty block structure. In this context, for = H* x (r],t)H y (rj,t), a complex matrix operator 
of input analytic IMFs to output analytic IMFs, pJ M,t) is a measure of the smallest uncertainty 
A = C 11X11 (note that if n v ^ n , then A = c' h ' XVx is not square) that causes instability if interpreted as 
a constant matrix feedback loop (ref. 28), 

^{a^^} (M) = o[H* x (Tl,t)HT(Ti,t)] 

and for scalar uncertainty structure A = { <5/7] : <5 e C} representing diagonal structure between 
input-output IMFs, 


n A (M) = p[H* x cn,t)H*cn,t)\ 


where p is the spectral radius (largest magnitude eigenvalue). This is distinctly different from 
the full-block structure in that uncertainty is only between corresponding input-output analytic 
IMFs (similar dominant frequencies), while ignoring uncertainty across different analytic IMFs 
(of different dominant frequencies) between input and output. This is obviously less conservative 
but generally less realistic as well, especially for nonlinear effects which cross frequencies. Larger 
values in either case represent effects of uncertainty between input and output such that higher 
correlated IMFs relate to more sensitivity to uncertainty at those dominant IMF frequencies. 

For MIMO signal analysis, it would be most appropriate to combine complex blocks (either 
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full or scalar) for each input-output into a multi-block structure, where each complex subblock 
corresponds to an input-output analytic IMF complex uncertainty structure. In computation this is 
often expanded in a block-diagonal context with repeated scalar and full blocks (ref. 28), where in 
this case each of these blocks would correspond to a single input-output analytic IMF uncertainty 
structure. 

To compare input-to-output correlations from aileron position to wing accels for different 
uncertainty structures, structured singular values were computed for the full-block and scalar-block 
structures between the collective aileron position and F-18 AAW wing accelerometer responses, 
then normalized with respect to the largest of the group. Table 4 lists the results for both types 
of uncertainty structures. Accelerometer#9 has the highest correlation ( = 1.00) with the input 
aileron position in either case. Interestingly, comparison between left and right wing accels is poor 
even from a symmetric aileron input, thereby indicating asymmetry and/or nonlinearity. Lower /U A 
values indicate a degree of robustness to uncertainty in a feedback stability sense, and as pointed 
out previously, scalar uncertainty is prevalently less conservative (as evidenced from absolute 
values before normalizing since the normalization factor is almost equal between the uncertainty 
structures). If an a priori bound on gain levels (normalization) between input-output is defined such 
that ||A|| x < 1 is considered acceptable in some sense (stability, health diagnostic, safety margin), 
then u A has an absolute (scaled) interpretation such that q A > 1 indicates that an acceptable 
threshold has been exceeded, and this can be used as a health, stability, or safety monitor. 


Table 4. Normalized results from F/A-18 active aeroelastic wing left- and right-wing 
acceleration data. 


Full 

3. L-wt-fwd 

7. L-wof-fwd 

11. L-wif-fwd 

9. R-wif-fwd 

5. R-wof-fwd 

1. R-wt-fwd 

Block 

4. L-wt-aft 

8. L-wof-aft 

12. L-wif-aft 

10. R-wif-aft 

6. R-wof-aft 

2. R-wt-aft 

Fwd-wing 

0.29 

0.29 

0.73 

1.00 

0.62 

0.80 

Aft-wing 

0.65 

0.67 

0.56 

0.19 

0.31 

0.72 


Scalar 

3. L-wt-fwd 

7. L-wof-fwd 

11. L-wif-fwd 

9. R-wif-fwd 

5. R-wof-fwd 

1. R-wt-fwd 

Block 

4. L-wt-aft 

8. L-wof-aft 

12. L-wif-aft 

10. R-wif-aft 

6. R-wof-aft 

2. R-wt-aft 

Fwd-wing 

0.11 

0.13 

0.19 

1.00 

0.47 

0.49 

Aft-wing 

0.23 

0.42 

0.27 

0.11 

0.06 

0.34 


A time-dependent interpretation is also available by calculating /r A (M,t) at each time point. 


[ * T 1 

H x ( rj, t t )H y ( r\, ti) V (scalar-block uncertainty structure) 
= (7 H*( Tj,ti )Hy (ri,tj ) \\/tj (full-block uncertainty structure). 


In figure 13 are plots of jU A M,t) for each input-to-output (aileron-to-accel) normalized analytic 
IMF where it is evident that values close to one (0 dB) are near a unity operator nonn limit. 
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Again, this depends on the bound A ^ <1 indicating an acceptable threshold, in this case maximum 
input-to-output gain normalized to one. Values close to one (0 dB) approach the acceptable limit. 
This uncertainty structure application has implications for model validation in the time domain 
(ref. 29). 


Aerostructures Test Wing 

Another example that includes an actual instability is taken from ATW (ref. 17) flight test data. 
The input is a sine sweep PZT voltage and outputs are three wingtip accelerometer responses. 
The EMDs of the PZT input and center wingtip accelerometer output near the flutter condition 
are displayed in figure 14. Flutter response EMDs of the center wingtip accelerometer and 
corresponding plots jU A ( M ,t ) between input analytic IMFs and all three wingtip accelerometer 
output analytic IMFs are shown in figure 15. The input PZT was not activated during the flutter 
occurrence, so the input is essentially a small arbitrary oscillation about zero. It is evident that at 
the point of instability past 7.5 s, (M,t) appropriately approaches a value of one (0 dB) since 
the input and output are normalized to unity, and assuming this limit corresponds to the limit of 
stability, and output response approaches its upper limit at flutter. 

PARAMETER ESTIMATION USING THE HILBERT-HUANG TRANSFORM 

The Hilbert transform can used for estimating modal parameters such as natural frequencies and 
damping ratios. For a single-mode system, the damped natural frequency can be easily determined 
from the Hilbert transform of the impulse response function. In this case, the damped natural 
frequency <x> d is given by the slope of the instantaneous phase angle plotted as a function of time, 
and damping ratio £ is estimated in a straightforward manner (ref. 30). 

For a linear single-mode system characterized by a pair of complex conjugate eigenvalues, 
the eigenvalues are directly related to the modal parameters as -£co n = <J d , where w n denotes the 
natural frequency. The impulse response function is given by 

h(t) = e ^ a>nt sin( co d t) 

niCQ d 


where m denotes the mass. Provided the damping ratio is small, the following relationships hold 
(ref. 30) 

H sin( co d t ) j = cos( (O d t ) 

H je -l =®" f cos( (O d t ) j = e _ ^®' if sin( co d t ) 
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where again H denotes the Hilbert transform. Therefore, the Hilbert transform of the impulse 
response function can be written as 

H{h(t)} := h(t) = — 1 — C os (co d t). 
m(O d 


The impulse response function can then be viewed as the real part of the following analytic signal 
(ref. 30) 


z(t) = h(t) + ih(t) 


= — *— e ^®" f sin (co d t) + i—^—e l =®' ,r 
m(O d m(O d 


cosf co d t ) 


The amplitude of this signal is given by 

a(t) = \]h(t) 2 +h( t ) 2 = — 1 — 

mco d 


which represents the envelope of the impulse response. Taking the natural logarithm of the 
amplitude yields 


In (a(t)) = -C,a> n t - ln( mco d ) 


The slope o of the natural log of a{t) versus time gives a = ~'C(o n . Given that co d can be obtained 
from the slope of the phase versus time, and using the definition of the damped natural frequency, 
there results co n , 




2 / 

=®n -cr 2 

= | 


Jo) d +(7 2 


and the damping ratio can be calculated as C,-o t co n . This approach to damping estimation is in 
the same spirit as the logarithmic decrement technique commonly used in structural dynamics. 
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Most systems have multiple modes, however, which implies that the impulse response function 
will generally contain contributions from several modes. Therefore, because they are valid only 
for single-mode systems, logarithmic decrement or the approach based on the Hilbert transform 
cannot be used for estimating damping. As emphasized earlier, the Hilbert transfonn does not 
yield meaningful frequency information for multiple-mode systems since it attempts to identify a 
single instantaneous frequency at each time step. The Hilbert-Huang algorithm makes it possible 
to apply the above technique for damping and frequency estimation because it can decompose a 
multi-component signal into a series of single-component signals through the EMD process. By 
performing EMD on the impulse response function of a multiple-mode system, the contributions 
of the different modes can be extracted and analyzed separately. This approach has been applied 
by Yang et al. (refs. 5-7) in a series of papers in which the Hilbert-Huang algorithm was used for 
parameter estimation of several multiple degree-of-freedom structures. A similar analysis is used 
here to estimate damping ratios and frequencies for aeroelastic systems. 


As an example, consider a prototypical linear pitch-plunge aeroelastic system that has been 
studied extensively in the literature. The system consists of an airfoil with pitch and plunge degrees 
of freedom, and the input is the deflection /3 of the trailing edge flap. The equations of motion for 
this system are given by 


m inx a b 
mx a b I a 



ch 0 
0 ch 



k h 0 


0 k, 


a 


= pU 2 b 


c l 

v n 


1 / 1 
( X H h H 

u u 


n 

A 

\ 

- a 

ba 

v 2 


y 


r 


i., in y 2 ' 


I h) 

\a I 


+ pU 2 b 


ba H — bh h — 


V 


u 


U 


— a 


V- 


b a 


-°Ir 


c h 


P 


where h denotes the plunge displacement of the airfoil and a represents the pitch angle. The terms 
on the left-hand side of the equation describe the structural dynamics of the system while the terms 
on the right-hand side represent quasi-steady aerodynamic forces and moments. In this study, the 
flow velocity U is varied until the linear flutter speed is reached at approximately U= 11.8 m/s. 
The values of the other system parameters are held constant as listed in table 5. 

Table 5. Parameters for the pitch-plunge aeroelastic system. 


m = 12.387 kg 

a = -0.6 

OC 

<N 

's 6 

ll 

k h = 2844.4 N/m 

I a = 0.065 m 2 kg 

x a = 0.2466 

c = -0.628 

m a 

k = 2.82 Nm 

a 

p =1.225 kg/m 3 

c h = 27.43 kg/s 

c, = 3.358 

l P 


b = 0.135 m 

c a = 0.036 m 2 kg/s 

c =-0.635 

m p 
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Because this system has two outputs, pitch and plunge, there are two impulse response functions 
that can be measured. Each of these functions generally contains contributions from both modes of 
the system. As an example, figure 16 depicts the pitch and plunge impulse response functions for 
a flow velocity of U = 8 m/s. These were obtained by simulating the pitch and plunge responses to 
an impulse applied to the trailing edge flap. The plunge impulse responses are used over a range 
of flow velocities to measure frequencies and damping ratios using the Hilbert-Huang algorithm. 
These impulse response functions are obtained by way of simulation in this report, but in practice 
they can be obtained by way of system identification techniques such as taking the inverse Fourier 
transform of the measured frequency response function or identifying first-order Volterra kernels 
from measured input-output data. 

The plunge impulse response functions were used to estimate the damped natural frequencies, 
damping ratios, and natural frequencies over a range of flow velocities. The results are summarized 
in table 6, which lists the actual and estimated damped natural frequencies and damping ratios for 
flow velocities ranging from U= 1 m/s to £/= 11.9 m/s (just past the linear flutter speed). The actual 
frequencies and damping ratios were calculated by solving for the eigenvalues of the system. The 
estimation of the modal parameters for U = 8 m/s is demonstrated in figures 17 and 18. Figure 17 
illustrates the EMD of the plunge impulse response function at U= 8 m/s. In this case, both modes 
of vibration are evident in the first IMF. This is because the first, higher-frequency mode damps 
out quickly and the second mode persists for a longer period of time. Therefore, the frequencies 
and damping ratios of both modes were estimated by using different portions of the first IMF, 
as depicted in figure 18. The figure shows plots of the IP and the natural log of the amplitude of 
the first IMF. Recall that the slopes of these lines are used to estimate the modal parameters. The 
initial portion of the amplitude plot is not a straight line because of initial transients in the IMF and 
possible boundary effects. Therefore, the portions of the plot that most resembled straight lines 
were used for the slope calculations for each mode, as labeled in the figure. Because this analysis 
technique is only approximate depending on the magnitude of the damping ratio (increased error 
for larger damping), the amplitude curve tends to oscillate about an average slope. Therefore, the 
slope in each case was obtained using a linear least-squares curve-fit of the data (refs. 5-7). 

The results presented in table 6 show that the estimates of the damped natural frequencies 
and damping ratios are fairly accurate in most cases. As the linear flutter speed (approximately 
U =11.8 m/s) is approached, the plunge mode is not discernible in the IMFs. In this regime, 
the system behaves essentially as a single-mode system and it was not possible to estimate the 
damping and frequency of the plunge mode. The actual and estimated damping ratios for each 
mode are plotted as a function of flow velocity in figure 19. 

A few observations are now presented. The first IMF was used to estimate the parameters of 
both modes. The analysis would be somewhat cleaner if each mode appeared in a different IMF. 
This is indeed the case in the work of Yang et al. (refs. 5-7) in which the frequency range of each 
IMF is controlled in the EMD process. Therefore, controlling the frequency range of each IMF is 
an option that could enhance the parameter estimation procedure. Also, even though parameters for 
both modes were estimated using a single IMF, the same results would generally not be obtained 
from an analysis of the impulse response function without the benefit of EMD. Although the IMF 
used in the analysis contains contributions from two modes, only one mode is present in the signal 
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at any point in time. In contrast, the impulse response function may contain contributions from 
both modes over the same time period. This example demonstrates some potential benefits of HHT 
for extracting modal parameters from aeroelastic flight data. 

Table 6. Actual and estimated damped frequencies and damping ratios. 


Velocity 

(m/s) 

Actual 

"/(Hz) 

Estimated 

"/Hz) 

Actual 

X, 

Estimated 

X 

1 

Plunge 

2.765 

2.735 

.0854 

.0874 


Pitch 

1.027 

1.022 

.0441 

.0444 

3 

Plunge 

2.739 

2.709 

.0855 

.0881 


Pitch 

1.057 

1.050 

.0503 

.0490 

5 

Plunge 

2.684 

2.655 

.0862 

.0895 


Pitch 

1.119 

1.096 

.0553 

.0591 

7 

Plunge 

2.596 

2.567 

.0877 

.0926 


Pitch 

1.217 

1.214 

.0590 

.0579 

9 

Plunge 

2.454 

2.452 

.0906 

.0957 


Pitch 

1.366 

1.362 

.0607 

.0637 

11 

Plunge 

2.198 

— 

.1012 

— 


Pitch 

1.628 

1.619 

.0544 

.0544 

11.5 

Plunge 

2.073 

— 

.1174 

— 


Pitch 

1.755 

1.744 

.0392 

.0394 

11.6 

Plunge 

2.040 

— 

.1262 

— 


Pitch 

1.788 

1.775 

.0305 

.0305 

11.7 

Plunge 

2.009 

— 

.1392 

— 


Pitch 

1.819 

1.806 

.0176 

.0017 

11.8 

Plunge 

1.984 


.1552 



Pitch 

1.844 

1.831 

.0016 

.0017 

11.9 

Plunge 

1.967 


.1713 



Pitch 


- 

-.0145 

.0146 
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CONCLUSIONS 


Application of the Hilbert-Huang algorithm for system signal decompositions, studying the 
effect of enhancements such as local-on-line behavior, understanding filtering properties, and 
especially for investigating correlations between input-output and between system sensors in 
terms of instantaneous properties, is revealed. System input-output signal analysis is introduced 
to characterize the time-varying amplitude and frequency components of multiple data channels, 
including input-to-output and distributed sensors, in terms of the intrinsic mode functions of 
the Hilbert-Huang transfonn. Significant departures from Fourier and other time-frequency or 
time-scale wavelet approaches are exposed. In these respects, this report attempts to show how the 
Hilbert-Huang transfonn behaves in a sometimes nonintuitive and subtle manner in the analysis 
of F-18 Active Aeroelastic Wing airplane aeroelastic flight test data. Online stability analyses and 
modal identification are also presented. Examples are given using aeroelastic test data from the F-18 
Active Aeroelastic Wing airplane, an Aerostructures Test Wing, and pitch-plunge simulation. 

An objective signal-adaptive basis function derivation, the Hilbert-Huang algorithm yields 
intrinsic mode functions giving instantaneous frequencies as functions of time that permit 
identification of imbedded structures. There is a multiresolution quality in the empirical mode 
decomposition process which even deals with intermittency by allowing multiple time-scales 
within an intrinsic mode function, but not allowing a similar time-scale simultaneously with other 
intrinsic mode functions. System identification in the intrinsic mode function subcomponent 
environment is a practical endeavor in the domain of multiresolution system identification. It 
should be noted that the idea of exploiting local properties for signal analysis applies to spatial data 
as well as temporal data with frequency and scale (translation and duration) variations. From the 
idea of empiquency to describe oscillations in images based on extrema points there are potential 
applications for general time-space-frequency-scale signal processing. 

Modern intelligent control and integrated aerostructures require control feedback signal 
processing cognizant of system stability and health. Time-varying linear or nonlinear modal 
characteristics derived from flight data are all within the realm of the Hilbert-Huang transform. 
Further research will investigate these issues and Hilbert-Huang transform connections between 
localized instantaneous dynamics, health diagnostics, and global system stability and performance 
for monitoring and prediction. 


Diyden Flight Research Center 

National Aeronautics and Space Administration 

Edwards, California, June 29, 2005 
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FIGURES 


Hilbert-Huang spectrum of x-| 



(a) Time-frequency spectra from the Hilbert transforms of the signal x y 
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(b) Time-frequency spectra from the Hilbert transform of the signal x 2 . 


Figure 1. Time-frequency spectra from the Hilbert transforms of the signals x t andx 2 . 
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Empirical mode decomposition of x-j 



(a) Empirical mode decomposition of the signal x T 


Hilbert-Huang spectrum of x-j 



(b) Hilbert-Huang spectrum of the signal x 2 . 


Figure 2. Empirical mode decomposition and the Hilbert-Huang spectrum of the signal x v 
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(a) Empirical mode decomposition of the signal x T 



(b) Hilbert-Huang spectrum of the signal x 2 . 


Figure 3. Empirical mode decomposition and the Hilbert-Huang spectrum of the signal x T 


31 


Empirical mode decomposition 






(a) Empirical mode decomposition. 


Signal Mode 1 
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Time 
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(b) Hilbert-Huang spectrum. 


Figure 4. Empirical mode decomposition and the Hilbert-Huang spectrum of the sum of two 
sinusoidal frequency modulations and one Gaussian wave packet. Original signal (top), FM1, 
FM2, Gaussian wave packet, and residual. 
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Empirical mode decomposition 
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Figure 5. Empirical mode decomposition of the sum of one sinusoid with two triangular 
waveforms. Original signal (top), high-frequency triangular, sinusoid, low-frequency triangular, 
and residual. 
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Empirical mode decomposition (constant plot scales) 





(a) Standard empirical mode decomposition. 


Empirical mode decomposition (constant plot scales) 




-IP? 






(b) Local-on-line empirical mode decomposition. 


Figure 6. Standard (a) and local-on- line (b) empirical mode decomposition of an F/A-18 active 
aeroelastic wing multisine aileron command input, with original signal at top and residual 
at bottom. 




Empirical mode decomposition for Accel#3 (left-fwd) 
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(a) Accel#3. 


Empirical mode decomposition for Accel#1 (right-fwd) 
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(b) Accel# 1. 


Figure 7. Empirical mode decompositions (with constant vertical plot scales) of F/A-18 active 
aeroelastic wing airplane left-fwd wingtip accel [(a), Accel#3], right-fwd wingtip accel [(b), 
Accel# 1], left-aft wingtip accel [(c), Accel#4], and right-aft wingtip accel [(d), Accel#2], responses 
from the multisine symmetric aileron command shown in figure 6 (original signal at top and 
residual at bottom in each set). 
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Empirical mode decomposition for Accel#4 (left aft) 






(c) Accel#4. 
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Empirical mode decomposition for Accel#2 (right aft) 
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(d) Accel#2. 


Figure 7. Concluded. 
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Figure 8. HLCC'l functions between intrinsic mode functions of Accel# 1 [figure 7(b)] and Accel#3 
[figure 7(a)], 
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(a) Instantaneous log-magnitude. 


Local transfer function phase 



200 

100 


lilij UjLj J Jill k.ikJlUllMliikla JJliLii j uliiliiLlililklilli.liJillili iuiiLUl .L.kii ilii . illLl Jl 1 Mi^Ui4ll1 AlJ 






M 


200 

100 

0 

200 

100 

0 

200 

100 


1 , "' i ' i '" 1, 1 1 1 W ; 




10 


15 


Time, s 


20 

050324 


(b) Instantaneous phase. 


Figure 9. Instantaneous log-magnitude [IM^, (a)] and phase [ IP ( X [L (b)] of intrinsic mode 
function transfer functions from F/A-18 active aeroelastic wing aileron input (figure 6) to Accel# 1 
[figure 7(b)]. 
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(a) Intrinsic mode, component-by-component. 



Empirigram 
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(b) Zoomed-in portion between 10 to 11 seconds. 


Figure 10. Empirical local correlation coefficient plots using two different representations of the 
empirical local correlation coefficient between the input signal from figure 6 to Accel# 1 in figure 
7(b). Top plot is an intrinsic mode component-by-component depiction, starting from the first 
at the top to the ninth at the bottom, where individual intrinsic functions are clearly delineated. 
Bottom plot is zoomed-in portion between 10 to 11 seconds of the same empirical local correlation 
coefficient. 


39 



Empirigram 
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Figure 11. Empirical local correlation coefficient represented using a three-dimensional contour of 
the result in figure 10. 
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Accel#3 


Accel#1 






(a) Instantaneous energy. 


Accel#3 


Accel#1 






(b) Marginal spectra. 

Figure 12. Instantaneous energy (a) and marginal spectra (b) of right-fwd wingtip accel (top 
right, Accel#l), right-aft wingtip accel (bottom right, Accel#2), left-fwd wingtip accel (top left, 
Accel#3), and left-aft wingtip accel (bottom left, Accel#4) responses, respectively in each set, 
from the multisine symmetric aileron command shown in figure 6. 
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(a) Time- varying for full-block uncertainy. 
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(b) Time-varying ,t) for scalar-blcok uncertainty. 


Figure 13. Time-varying j,i A ( M,t) (dB) for full-block (a) and scalar-block (b) uncertainty. 
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Empirical mode decomposition (constant plot scales) 
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(a) Empirical mode decomposition of the aerostructures test wing input PZT. 


Empirical mode decomposition (constant plot scales) 
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(b) Empirical mode decomposition of the center wingtip accelerometer. 


Figure 14. Empirical mode decomposition of the aerostructures test wing input PZT (a) and center 
wingtip accelerometer (b) at Mach 0.82 and an altitude of 10,000 ft (3,048 m) just before flutter 
(original signal at top and residual at bottom in each set). 
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Empirical mode decomposition (constant plot scales) 
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(a) Empirical mode decomposition of the center wingtip accelerometer. 
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(b) Time-varying for full-block uncertainty between input and output intrinsic mode 

functions of the three wingtip accelerometer outputs. 

Figure 15. Empirical mode decomposition of the center wingtip accelerometer (a) near Mach 0.83 
and an altitude of 10,000 ft (3,048 m) at flutter (original signal at top and residual at bottom), 
and (b) time-varying ji^(M,t) for full-block uncertainty between input and output intrinsic mode 
functions of the three wingtip accelerometer outputs. 
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x 10 -4 Impulse response, h - Plunge, U = 8 m/s 



(a) Plunge and pitch; impulse response, h. 



(b) Plunge and pitch; impulse response, a. 


Figure 16. Plunge and pitch impulse response functions for U = 8 m/s. 
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Figure 17. Empirical mode decomposition of the plunge impulse response function 
(m/rad) at U = 8m/s. 
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(a) Phase angle. 



(b) Amplitude. 


Figure 18. Instantaneous phase and amplitude plots used to estimate modal parameters from the 
first intrinsic mode function of the plunge impulse response function at 17 = 8 m/s. 
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Figure 19. Actual and estimated damping ratios. 
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